#include <iostream>
#include <fstream>
#include <string>
#include "fields.h"
#include "cellinf.h"
#include "simbox.h"
#include "laser.h"
#include "common.h"
#ifdef __MPI
#include "mpi_data.h"
#endif
#ifdef __OMP
#include "omp.h"
#endif
#include "injector.h"
#include "denfunc.h"
#include "detector.h"
using namespace std;

long __Ptcls::total_id = 0;
long __Ptcls::total_num = 0;
long& __Injector::pid = __Ptcls::total_id;
int main(int argc, char* argv[])
{

    int myrank = 0, ranksize = 1;
#ifdef __MPI
    MPI_Status status;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    MPI_Comm_size(MPI_COMM_WORLD, &ranksize);
#endif

#ifdef __OMP
    omp_set_num_threads(4);
#endif
    /***
    __Ptcls myptcls;
    cout<<"size of ptcls is "<<sizeof(myptcls)<<endl;
    return 0;
    ***/
    string simname = simu_name;
	__Date simdate;
    int x_cell_sim, y_cell_sim;
    double x_len_sim, y_len_sim;
    x_cell_sim = x_cell; y_cell_sim = y_cell;
    x_len_sim= x_len; y_len_sim = y_len;
    __Vect2<double> time_range;
    double dx = x_len_sim / x_cell_sim;
    double dy = y_len_sim / y_cell_sim;
    double delta_x_box = x_len_sim / ranksize, delta_y_box = y_len_sim;
    __Vect2<double> dxy(dx, dy);


    // 有两个是ghost cell
    double xmin_box = delta_x_box * myrank - dx, xmax_box = delta_x_box * (myrank + 1) + dx;
    double ymin_box = - dy, ymax_box = delta_y_box + dy;
    
    // ghost cell size

    __Vect2<double> basexymax(x_len_sim, y_len_sim);
    __Vect2<double> basexymin(0, 0);
    
	
    __Vect2<int> xy_cell_num(x_cell_sim / ranksize  + 2, y_cell_sim + 2);
    // ghost cell
    
    int x_grid_min = xy_cell_num.member[0] * myrank - myrank, x_grid_max = \
                      xy_cell_num.member[0] * (myrank + 1) - myrank;
    int y_grid_min = 0, y_grid_max = \
                      xy_cell_num.member[1];
    __Vect2<int> xygridmin(x_grid_min, y_grid_min);
    __Vect2<int> xygridmax(x_grid_max, y_grid_max);
	
    //cout<<myrank<<" grid min "<<xygridmin<<" max "<<xygridmax<<endl;
    //cout<<" cell num "<<xy_cell_num<<endl;
    __Vect2<double> xrange(xmin_box, xmax_box);
	__Vect2<double> yrange(ymin_box, ymax_box);
	
	__Vect2<double> xymin(xmin_box, ymin_box);
    __Vect2<double> xymax(xmax_box, ymax_box);
    
    string fname = file_name;
    __Vect2<double> pos(dxy.member[0], y_len_sim * 0.5); 
    double gamma = 20; // 10 MeV
    double vx = c_light_speed * sqrt(1.0 - 1.0 / gamma / gamma);
    __Vect3<double> velocity(vx, 0, 0.0);
    /// species to be merged
    vector<int> msp = MSP;
    int max_ppc = MAX_PPC;
    double delta_e = MERGE_BIAS_E;
    double delta_v = MERGE_BIAS_V;
    //cout<<"rank size = "<<ranksize<<endl;
    __Vect2<int> position(myrank, 0);
    cout<<myrank<<" xy min "<<xymin<<"\t"<<" xy max "<<xymax<<" dxy "<<dxy<<endl;
    
    // 没有设置，以后在设置
    int bd_type[4] = bc_type;
    /***
    if(rankyj == 0)
    {
        bd_type[0] = 1;
    }
    if(rankyj == 1)
    {
        bd_type[1] = 1;
    }
    ***/
    cout<<myrank<<"bd typ "<<bd_type[0]<<"\t"<<bd_type[1]<<"\t"<<\
        bd_type[2]<<"\t"<<bd_type[3]<<endl; 
    vector<string> fname1d, fname2d;
    string sigmabr, sigmabh;
    sigmabr = SIGMABR;
    sigmabh = SIGMABH;
    fname1d.push_back("tables/etaforp.txt");
    fname1d.push_back("tables/chiforp.txt");
    fname1d.push_back("tables/fforp.txt");
    fname1d.push_back("tables/heta.txt");
    fname1d.push_back("tables/tchi.txt");
    fname1d.push_back("tables/chifort.txt");
    fname1d.push_back("tables/etaforh.txt");
    fname1d.push_back("tables/ele_ek.txt");
    fname1d.push_back("tables/xrange.txt");
    fname2d.push_back("tables/petachi.txt");
    fname2d.push_back("tables/pchif.txt");
    fname2d.push_back("tables/chietap.txt");
    fname2d.push_back(sigmabr);
    fname2d.push_back(sigmabh);
    int startind = damp_length;
    // 下上左右
    int startindex[4] = {startind, int(y_cell_sim) - startind, startind, \
        int(x_cell_sim) / ranksize - startind};
    //cout<<"startindex "<<startindex[0]<<"\t"<<startindex[1]<<"\t"\
        <<startindex[2]<<"\t"<<startindex[3]<<endl;
    int damp_len[4] = {startind, startind, startind, startind};
	int species = species_all;
	double delta_t;
    delta_t = dx / c_light_speed / sqrt(2.0) * 0.8;
    if(DISABLE_EM)
    {
#ifdef DELTA_T
        delta_t = DELTA_T;
#endif
    }
    int t_steps = (time.member[1] - time.member[0]) / delta_t;
    double wavelength = wave_length;
    double duration = duration_um; // cTl = 3.2 um
    double t_center = laser_center;
    double waist = laser_waist;
    double scalea = scaled_a;
    double omega = laser_omega;
    double intensity = laser_intensity;
    double amp = laser_amplitude;
    int polarization = polarize; // z polarized
    double dens = DENSITY;
    double mydens = dens;
    double dens1[4] = {mydens, mydens / ATOM_NUM, 0, 0};
    double dens3[4] = {0, 0, 0, 0};
    double dumpsec = dump_in_sec;
    int dumper = dumpsec / delta_t;
    int ifdump = DUMP;
    //double amp = get_intensity(20, omega); 
    if(myrank == 0)
    {
        /***
        cout<<" density "<<dens<<endl;
        cout<<" a = "<<scalea<<" so, Intensity = "<<intensity<<" w/m^2"<<endl;
        cout<<" amp = "<<amp<<endl;
        cout<<" v dt "<<velocity.member[0] * delta_t<<endl;
        cout<<" dx "<<dx<<endl;
        cout<<" corant condition "<<(dx > c_light_speed * delta_t)<<endl;
        cout<<" steps: "<<t_steps<<endl;
        cout<<" dx / v = "<<dx / velocity.member[0] / delta_t<<endl;
        cout<<" wavelength = "<<wavelength<<endl;
        cout<<" dens of cell "<<mydens<<endl;
        cout<<" period = "<<wavelength / dx<<endl;
        cout<<" need dt "<<wavelength / c_light_speed / delta_t<<endl;
        cout<<" dx / c / dt = "<<dx / c_light_speed / delta_t<<endl;
        cout<<" omegape * delta_t = "<<omegape(dens) * delta_t<<endl;
        cout<<" delta_t = "<<delta_t<<endl;
        cout<<" xy cell num = "<<xy_cell_num<<endl;
        // 两个约束条件，corant condition，omegape * dt < 2
        if(dx < c_light_speed * delta_t)
        {
            return 1;
        }
        ***/
    }
    cout<<myrank<<" dumper = "<<dumper<<" delta t "<<delta_t<<endl;
    /// 0 不需要开启，１需要开启, 顺序0-syn, 1-brem, 2-bwpair, 3-bhpair
    int monte_flag[4] = mc_flag;
	int maxid = 0;
    vector<Particles> spptcls;
    vector<int> ppc = {};
#ifdef __INITSP
    init_sp(spptcls, ppc);
#endif
    cout<<"ppc size "<<ppc.size()<<endl;
    //cout<<"rank "<<myrank<<" will init box"<<endl;
    __Ptcls tempptcls;
    int local_init_id = max_ppc * species * xy_cell_num.member[0] * xy_cell_num.member[1] * myrank; 
    tempptcls.total_id = local_init_id;
    __Simbox mybox(simname, simdate, \
			xygridmin, xygridmax, xy_cell_num, xrange, \
			yrange, time,t_steps, delta_t, species, spptcls,\
            maxid, ppc, position, DISABLE_EM);
    mybox.set_data_file(fname1d, fname2d);
    mybox.set_monte_carlo(monte_flag);
    mybox.set_field_boundary(bd_type, startindex, damp_len); 
    mybox.init();
    cout<<"box init over "<<endl;
    mybox.max_ppc = max_ppc;
    mybox.delta_e = delta_e;
    mybox.delta_v = delta_v;
    __Cellgroup& mygroup = mybox.cellgroup;
    //density_dist();
    //cout<<"myrank is "<<myrank<<" myrank "<<myrank<<" per "<<ranksize<<endl;
#ifndef __USRDENS
    if(myrank == 0)
    {

        for(int i = 1; i < mygroup.member.size() - 1; i ++)
        {
            for(int j = 1; j < mygroup.member[i].size() - 1; j ++)
            {
                //if(i > 0 && i <32)
                if(i > 10)
                {
                    mygroup.member[i][j].set_density(dens1);
                }
            }
        }
    }
    else if(myrank == ranksize - 1)
    {
        cout<<"myrank is - "<<myrank<<"\t"<<endl;
        
        for(int i = 1; i < mygroup.member.size() /2 ; i ++)
        {
            for(int j = 1; j < mygroup.member[i].size() - 1; j ++)
            {
                mygroup.member[i][j].set_density(dens1);
            }
        }
        for(int i = mygroup.member.size() / 2; i < mygroup.member.size() ; i ++)
        {
            for(int j = 1; j < mygroup.member[i].size() - 1; j ++)
            {
                mygroup.member[i][j].set_density(dens3);
            }
        }
    }
    else
    {
        for(int i = 1; i < mygroup.member.size() - 1; i ++)
        {
            for(int j = 1; j < mygroup.member[i].size() - 1; j ++)
            {
                mygroup.member[i][j].set_density(dens1);
            }
        }
    }
#endif
    cout<<myrank<<" initing ptcls"<<endl; 
    mygroup.init_ptcls();
    cout<<" ptcls init over "<<endl;
    //mygroup.set_sp_velocity(0, velocity); 
    mybox.init_flag();
    //cout<<myrank<<"'s min is "<<mygroup.member[0][0].xymin<<" max is "\
        <<mygroup.member[mygroup.member.size() - 1][mygroup.member[0].size() - 1].xymax<<endl;
    cout<<" group size "<<mygroup.member.size()<<"\t"<<mygroup.member[0].size()<<endl;
    cout<<" field size "<<mybox.fieldgroup.member.size()<<"\t"<<mybox.fieldgroup.member[0].size()<<endl;
    __Vect2<int> gridsize(mybox.fieldgroup.member.size(), mybox.fieldgroup.member[0].size());
    __Laser mylaser(0, 0, delta_t, t_center, duration, pos, wavelength, xymin, xymax,\
                gridsize, dxy, waist, amp, polarization);
    cout<<"laser has been created "<<endl;
    __Vect2<int> ipos(1, y_cell_sim / 2);
    long total_inject = 1e7;
    cout<<"ipos "<<ipos<<"\t"<<y_len<<endl;
    __Injector mygun1(0, 0, -e_charge, e_mass, ipos, velocity, total_inject, \
            1.0e-6, DENSITY, dx, dy);    
    __Vect2<double> gps1(dx, y_cell_sim / 2 * dy);
    mygun1.gps = gps1;

    __Vect2<double> ekrange(MeV, 10.0 * MeV);
    __Vect2<double> detect_xymin(100 * UNITS, 2 * UNITS);
    __Vect2<double> detect_xymax(101 * UNITS, 150 * UNITS);
    __Detector mydetector(0, 3, ekrange, detect_xymin, detect_xymax, dxy);
    mydetector.init(mybox);


    /// open a file for time step out put
    string tfile = fname + "time.txt";
    ofstream myout;
    myout.open(tfile);
    for(int i = 0; i < t_steps; i++)
    {
        long pdata_tags[8] = {i, t_steps + i, 2 * t_steps + i, 3 * t_steps + i, \
            4 * t_steps + i, 5 * t_steps + i, 6 * t_steps + i, 7 * t_steps + i};
        long fdata_tags[16] = {\
            8 * t_steps + i, 9 * t_steps + i, 10 * t_steps + i, 11 * t_steps + i, \
            12 * t_steps + i, 13 * t_steps + i, 14 * t_steps + i, 15 * t_steps + i, \
            16 * t_steps + i, 17 * t_steps + i, 18 * t_steps + i, 19 * t_steps + i, \
            20 * t_steps + i, 21 * t_steps + i, 22 * t_steps + i, 23 * t_steps + i};
        // 向左发送6个，三个size，三个data
        // 向右同，先Field，后ptcls, 再cellfield
        
        if(myrank == 0)
        {
            cout<<" step : "<<i<<endl;
            //cout<<" ptcls id max "<<tempptcls.total_id<<"\t"<<tempptcls.total_num<<endl;
        }
        /***
        if(mygun1.gps < mybox.xymax && mybox.xymin < mygun1.gps)
        {
            mygun1.trigger(mybox);
        }
        ***/
        if(myrank == 0 && DISABLE_EM == 0)
        {
            mylaser.update();
            mylaser.setting_fields(mybox);
        }
        /// 第一步是更新电磁场的B


#ifdef __MPI
        mpi_exchange_fields_lr(myrank, ranksize, mybox, fdata_tags, status);
        MPI_Barrier(MPI_COMM_WORLD); 
        shift_data_tags(fdata_tags);
#endif
#ifdef __Debug
        cout<<myrank<<" before advance 1 "<<endl;
#endif
        mybox.advancest1(myrank);
#ifdef __Debug
        cout<<myrank<<" after advance 1 "<<endl;
#endif
        /// 交换数据,保证相邻盒子的B也被更新
#ifdef __MPI
        /// 更新电场
        mpi_exchange_fields_lr(myrank, ranksize, mybox, fdata_tags, status);
        MPI_Barrier(MPI_COMM_WORLD); 
        shift_data_tags(fdata_tags);
#endif
#ifdef __Debug
        cout<<myrank<<" before advance 2 "<<endl;
#endif
        mybox.advancest2(myrank);
#ifdef __Debug
        cout<<myrank<<" after advance 2 "<<endl;
#endif
        /// 更新磁场后半步, 将电磁场拷贝到cell中，并weight到粒子上
#ifdef __MPI 
        mpi_exchange_fields_lr(myrank, ranksize, mybox, fdata_tags, status);
        MPI_Barrier(MPI_COMM_WORLD); 
        shift_data_tags(fdata_tags);
#endif
#ifdef __Debug
        cout<<myrank<<" before advance 3 "<<endl;
#endif
        mybox.advancest3(myrank);
#ifdef __Debug
        cout<<myrank<<" after advance 3 "<<endl;
#endif

        /// 交换数据，保证电磁场更新完全,粒子可以利用这部分场运动
        /// 更新粒子，运动，会出现越界粒子
#ifdef __MPI   
        mpi_exchange_fields_lr(myrank, ranksize, mybox, fdata_tags, status);
        MPI_Barrier(MPI_COMM_WORLD); 
        shift_data_tags(fdata_tags);
#endif
#ifdef __Debug
        cout<<myrank<<" before advance 4 "<<endl;
#endif
        mybox.advancest4(myrank, i);
#ifdef __Debug
        cout<<myrank<<" after advance 4 "<<endl;
#endif
        /// 交换数据，保证越界粒子被发送到隔壁盒子的越界粒子束中
#ifdef __MPI
        mpi_exchange_fields_lr(myrank, ranksize, mybox, fdata_tags, status);
        MPI_Barrier(MPI_COMM_WORLD); 
        shift_data_tags(fdata_tags);
        //cout<<myrank<<" begin ex ptcls"<<endl;
        
        mpi_exchange_ptcls(myrank, ranksize, mybox, pdata_tags, status);
        MPI_Barrier(MPI_COMM_WORLD); 
#endif
        /// 将越界粒子发送到自己的网格，将粒子产生的电荷与电流weight到网格上,然后拷贝到field中
#ifdef __Debug
        cout<<myrank<<" before advance 5 "<<endl;
#endif
        mybox.advancest5(myrank);
#ifdef __Debug
        cout<<myrank<<" after advance 5 "<<endl;
#endif
        //mydetector.detect(mybox); 
        /// 交换数据，保证网格,场数据交换完成
        //cout<<myrank<<" over detect cell copy to field "<<endl;
        mybox.cell_copy_to_field(0, myrank);
        //cout<<"copy over"<<endl;
        if(i % 10 == 0)
        {
            for(int imsp = 0; imsp < species; imsp ++)
            {
                if(msp[imsp] == 1)
                {
                    //cout<<"begin merging "<<imsp<<endl;
                    mybox.merge(imsp);
                    //mybox.wide_merge(imsp, 2);
                    //cout<<"end of merging"<<endl;
                }
            }
        }
        if(i%dumper == 1 && ifdump == 1)
        {
            //cout<<"myrank "<<myrank<<" maxid "<<tempptcls.total_id<<endl;
            if(myrank == 0)
            {
                myout<<i<<"\t"<<delta_t * i<<endl;
            }
            cout<<" "<<myrank<<" dumping "<<i<<endl;
            if(!DISABLE_EM)
            {
                mybox.fieldgroup.dump(myrank, i, fname, 0);
                // 0 rho, 1 ex, 6 bz, 7 jx
                mybox.fieldgroup.dump(myrank, i, fname, 1);
                mybox.fieldgroup.dump(myrank, i, fname, 2);
                mybox.fieldgroup.dump(myrank, i, fname, 3);
                mybox.fieldgroup.dump(myrank, i, fname, 4);
                mybox.fieldgroup.dump(myrank, i, fname, 5);
                mybox.fieldgroup.dump(myrank, i, fname, 6);
                mybox.fieldgroup.dump(myrank, i, fname, 7);
                mybox.fieldgroup.dump(myrank, i, fname, 8);
                mybox.fieldgroup.dump(myrank, i, fname, 9);
                // 0-e- 1 ion 2 e+ 3 gamma
                
                mybox.cellgroup.dump(myrank, i, fname, 0);
                mybox.cellgroup.dump(myrank, i, fname, 1);
                mybox.cellgroup.dump(myrank, i, fname, 2);
                mybox.cellgroup.dump(myrank, i, fname, 3);
                // ek e- 4, ion 5, pos 6
                //mybox.cellgroup.dump(myrank, i, fname, 4);
            
            }
            mybox.cellgroup.dump(myrank, i, fname, 5);
            mybox.cellgroup.dump(myrank, i, fname, 6);
            mybox.cellgroup.dump(myrank, i, fname, 7);
            mybox.cellgroup.dump(myrank, i, fname, 8);
            //mybox.cellgroup.dump(myrank, i, fname, 9);

            mybox.cellgroup.dump_ptcls(myrank, i, fname, 0);
            mybox.cellgroup.dump_ptcls(myrank, i, fname, 1);
            mybox.cellgroup.dump_ptcls(myrank, i, fname, 2);
            mybox.cellgroup.dump_ptcls(myrank, i, fname, 3);
            //mybox.cellgroup.dump_ptcls(myrank, i, fname, 4);
            //mydetector.dump(myrank, i, fname);
            //mydetector.reset();

        }
        if(i == EXIT)
        {
            //mydetector.dump(myrank, i, fname);
            myout.close();
#ifdef __MPI
            MPI_Finalize();
#endif
            return 0;
        }
    }
    myout.close();
    //mydetector.dump(myrank, t_steps, fname);
#ifdef __MPI
    MPI_Finalize();
#endif
    return 0;
}


